# Descriptive stats 
skim(df)

plot_ly(data = df, 
        x = ~profit,
        type = "histogram",
        histnorm = "probability") %>%  
  layout(title = "Frequency distribution of Profit",
         xaxis = list(title = "Profit",
                      zeroline = FALSE),
         yaxis = list(title = "Frequency",
                      zeroline = FALSE))

# Response variable renaming  
df_model_hot <- df_model_hot %>%
  rename(Y = profit)

# Separate into random training and test data sets ----

# set.seed() for reproducibility 
set.seed(42) 

# We keep 80% of the data for model training 
split_strat <- initial_split(df_model_hot, prop = 0.8)
df_train  <- training(split_strat)
df_test   <- testing(split_strat)

# Training set 
X <- as.data.table(df_train %>% select(-Y))
Y <- df_train$Y

XY <- cbind(Y,X)
names(XY)[1] <- "Y"

# Testing set 
X_test <- as.data.table(df_test %>% select(-Y))
Y_test <- df_test$Y

XY_test <- cbind(Y_test, X_test)
names(XY_test)[1] <- "Y"

# Cross validation 
# Higher CV metrics determine that our model does not suffer from high variance and generalizes well on unseen data

# Number of folds
N <- 5

# Folds building
set.seed(42)
Folds <- sample(1:dim(XY)[1] %% N + 1)

Output <- list()
Time <- list()

i <- 1

# for loop to perform cross validation 
for(i in 1:N){
  Predictions <- list()
  
  # Estimation samples
  XE <- X[Folds != i]
  YE <- Y[Folds != i]
  # Data needs to be in a matrix format for some models 
  XME <- as.matrix(XE)
  YME <- as.matrix(YE)
  XYE <- XY[Folds != i]
  
  # Prediction samples 
  XP <- X[Folds == i]
  YP <- Y[Folds == i]
  # Data needs to be in a matrix format for some models 
  XMP <- as.matrix(XP)
  YMP <- as.matrix(YP)
  XYP <- XY[Folds == i]
  
  print(stringr::str_c("Fold number ", i))
  print("**************************************************")
  
  # Add the Y variable to the prediction dataset (for comparison purposes)
  Predictions[["Y"]] <- YP
  Predictions[["Y"]] <- as.numeric(Predictions[["Y"]])
  
  #********************* ORDINARY LEAST SQUARES **********************#
  
  print(Nom <- "OLS (LM)")
  Start <- proc.time()
  Model <- stats::lm(formula = Y ~ ., data = XYE)
  Response <- predict(Model, newdata = XP, type = "response")
  Predictions[[Nom]] <- as.numeric(Response)
  Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  
  #********************* REGULARIZED REGRESSIONS *********************#
  
  print(Nom <- "OLS (GLMNET)")
  Start <- proc.time()
  Model <- glmnet::glmnet(x = XME, y = YME, lambda = 0, alpha = 0, family = "gaussian")
  Response <- predict(Model, newx = XMP, type = "response")
  Predictions[[Nom]] <- as.numeric(Response)
  Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  
  print(Nom <- "Ridge (GLMNET)")
  Start <- proc.time()
  Model <- glmnet::cv.glmnet(x = XME, y = YME, alpha = 0, family = "gaussian")
  Response <- predict(Model, s = "lambda.min", newx = XMP, type = "response")
  Predictions[[Nom]] <- as.numeric(Response)
  Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  
  # Regularized Regression to reduce the model variance 
  # alpha = 1 for lasso (push some coefficients to 0 for some variables)
  # alpha = 0 for ridge (no coefficients at 0)
  
  # λ is a tuning parameter that helps to control our model from over-fitting to the training data
  # However, to identify the optimal λ value we need to perform cross-validation (CV)
  # cv.glmnet provides a built-in option to perform k-fold CV, and by default, performs 10-fold CV
  
  # the 10-fold CV mean squared error (MSE) across the λ values
  
  # number of the top = number of variables 
  # The first vertical dashed lines represent the λ value with the minimum MSE 
  # and the second one the largest λ value within one standard error of the minimum MSE
  
  # min(Model$cvm)       # minimum MSE
  # Model$lambda.min     # lambda for this min MSE
  
  # Model$cvm[Model$lambda == Model$lambda.1se]  # 1 st.error of min MSE
  
  #  Model$lambda.1se  # lambda for this MSE
  
  # Influential variables 
  
  # coef(Model, s = "lambda.1se") %>%
  #   tidy() %>%
  #   filter(row != "(Intercept)") %>%
  #   ggplot(aes(value, reorder(row, value), color = value > 0)) +
  #   geom_point(show.legend = FALSE) +
  #   ggtitle("Influential variables") +
  #   xlab("Coefficient") +
  #   ylab(NULL)
  
  print(Nom <- "Lasso (GLMNET)")
  Start <- proc.time()
  Model <- glmnet::cv.glmnet(x = XME, y = YME, alpha = 1, family = "gaussian")
  Response <- predict(Model, s = "lambda.1se", newx = XMP, type = "response")
  Predictions[[Nom]] <- as.numeric(Response)
  Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  
  # coef(Model, s=Model$lambda.1se)
  
  Alphas <- c(0.25, 0.5, 0.75)
  
  for(Alpha in Alphas){
    print(Nom <- stringr::str_c("Elastic net (GLMNET|alpha=", Alpha, ")"))
    Start <- proc.time()
    Model <- glmnet::cv.glmnet(x = XME, y = YME, alpha = Alpha, family = "gaussian")
    Response <- predict(Model, s = "lambda.min", newx = XMP, type = "response")
    Predictions[[Nom]] <- as.numeric(Response)
    Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  }
  
  #********************* TREE *******************************#
  
  print(Nom <- "Tree (RPART)")
  Start <- proc.time()
  Model <- rpart::rpart(formula = Y ~ ., data = XYE)
  Response <- predict(Model, newdata = XP, type = "vector")
  Predictions[[Nom]] <- as.numeric(Response)
  Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  
  # rpart.plot(Model)
  
  #********************* RANDOM FORESTS *********************#
  # tune the number of tree hyperparameter in the cross validation 
  Numbers_Of_Trees <- c( 100, 250, 500)
  
  # make ranger compatible names
  names(XYE) <- make.names(names(XYE), allow = FALSE)
  names(XP) <- make.names(names(XP), allow = FALSE)
  
  for (ntree in Numbers_Of_Trees) {
    
    print(Nom <- stringr::str_c("Random forest (RANDOMFOREST|ntree=", ntree, ")"))
    Start <- proc.time()
    features <- setdiff(names(XYE), "Y")
    Model <- ranger(formula    = Y ~ .,
                    data       = XYE,
                    num.trees  = ntree, # number of trees 
                    mtry       = floor(length(features) / 3), # the number of variables 
                    # to randomly sample as candidates at each split
                    # usal value for regression = p/3 with p the total number of variables  
                    respect.unordered.factors = 'order',
                    verbose    = FALSE,
                    seed       = 42)
    Response <- predict(Model, data = XP, type = "response")
    Response <- Response$predictions
    Predictions[[Nom]] <- as.numeric(Response)
    Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  }
  
  #********************* BOOSTING *************************#
  
  # Xgboost - apprx 10x faster than gbm
  # For boosting tree models, the most important parameter is the number of trees to build - nrounds parameter
  # XGBoost only works with matrices that contain all numeric variables
  
  # tune the number of tree hyperparameter in the cross validation 
  Numbers_Of_Trees <-  c(100, 250, 500)
  
  for (nrounds in Numbers_Of_Trees){
    print(Nom <- stringr::str_c("XGBoost|n.trees=", nrounds))
    Model <- xgboost(data = XME, 
                     label = YME,
                     nrounds = nrounds,
                     eta = 0.1, # learning rate 0.3 by default
                     max_depth = 4, #  Maximum depth of a tree. 6 by default
                     objective = "reg:linear",
                     print_every = 1000) # hide intermediate results 
    Response <- predict(Model , XMP, missing = NA)
    Predictions[[Nom]] <- as.numeric(Response)
    Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  }
  
  # #********************* SVM *****************************#
  # NOT RUN 
  # tune Kernel and Cost hyper-parameters in the cross validation 
  # Kernels <- c("linear", "polynomial", "radial")
  # Costs <- c(0.1, 1, 2)
  # 
  # for(Kernel in Kernels)
  # {
  #   for(Cost in Costs)
  #   {
  #     print(Nom <- stringr::str_c("SVM (e1071|kernel=", Kernel, "|cost=", Cost, ")"))
  #     Start <- proc.time()
  #     Model <- e1071::svm(Y ~ ., data = XYE, kernel = Kernel, cost = Cost)
  #     Response <- predict(Model, newdata = XP, type = "response")
  #     Predictions[[Nom]] <- as.numeric(Response)
  #     Time[[Nom]] <- ifelse(is.null(Time[[Nom]]), 0, Time[[Nom]]) + (proc.time() - Start)[3]
  #   }
  # }
  
  Output[[1 + length(Output)]] <- as.data.table(do.call(cbind, Predictions))
  
  print("**************************************************")
}

# Stack all the outputs tables up 
Output <- data.table::rbindlist(l = Output, use.names = TRUE, fill = TRUE)

# Measure execution time 
Time <- data.table(Model = names(Time), Time = unlist(Time))

# Function to define the error metric : Mean Absolute Error 
loss_function <- function(X,Y){
  mean(abs(X-Y))
}

# Assess cross validation errors for each model 
Errors <- apply(X = Output, MARGIN = 2, FUN = loss_function, Y = Output$Y)
Errors <- data.table(Model = names(Errors), error_rate = Errors)
Errors <- merge(x = Errors, y = Time, by = "Model", all.x = TRUE, all.y = FALSE)
Errors <- Errors[order(error_rate)]
Errors

# Training set 
X <- as.data.table(df_train %>% select(-Y))
Y <- df_train$Y

XY <- cbind(Y,X)
names(XY)[1] <- "Y"

# Testing set 
X_test <- as.data.table(df_test %>% select(-Y))
Y_test <- df_test$Y

XY_test <- cbind(Y_test, X_test)
names(XY_test)[1] <- "Y"

# make ranger compatible names
names(XY) <- make.names(names(XY), allow = FALSE)
names(X_test) <- make.names(names(X_test), allow = FALSE)

features <- setdiff(names(XY), "Y")

# Best Model: Random Forest 500 trees according to the the cross validation errors 

Model_final<- ranger(formula    = Y ~ .,
                     data       = XY,
                     num.trees  = 500, 
                     mtry       = floor(length(features) / 3), 
                     respect.unordered.factors = 'order',
                     # To plot variable imortance 
                     importance = "impurity",
                     # Quantile regressions forest to assess prediction intervals
                     quantreg   = TRUE,
                     verbose    = FALSE,
                     seed       = 42)

# Predict on the test set
predictions_testRF <- predict(Model_final,
                              data = X_test, 
                              type = "response")

predictions_test <- data.frame(cbind(Y_test, predictions_testRF$predictions))

# Performance on the test set
loss_function(X = predictions_testRF$predictions, Y = Y_test)

#  Plot Top important variables 

var_imp_ranger<-as.data.frame(Model_final$variable.importance)
variables <-(as.vector((row.names(var_imp_ranger))))
var_imp_ranger <- cbind(variables,var_imp_ranger)
var_imp_ranger <- var_imp_ranger %>% 
  arrange(desc(`Model_final$variable.importance`))


p <- ggplot(var_imp_ranger,
       aes(x = reorder(variables,`Model_final$variable.importance`) ,
           y = `Model_final$variable.importance`, 
           fill = `Model_final$variable.importance`))+ 
  geom_bar(stat="identity",
           position="dodge",
           aes(text = `Model_final$variable.importance`))+
  coord_flip()+
  ylab("Importance")+
  xlab("")+
  ggtitle("Variable Importance - Random Forest 500 Trees")+
  guides(fill = F)+
  scale_fill_gradient()+
  theme_classic()

ggplotly(p, tooltip = "text")

# Quantile Prediction
quantiles_reg_forests <- predict(Model_final, 
                                 data = X_test,
                                 type = "quantiles",
                                 quantiles = seq(0, 1, by = 0.1))

pred_intervals <- as.data.table(quantiles_reg_forests$predictions)

# Concatenate actual value (profit) and predicted quantiles  
actual_profit_w_pred_intervals <- data.table(cbind(Y_test, pred_intervals))

ncols <- ncol(actual_profit_w_pred_intervals)

names(actual_profit_w_pred_intervals)[2:ncols] <- c("Q0", "Q10","Q20", "Q30", "Q40", "Q50",
                                                    "Q60", "Q70","Q80", "Q90", "Q100")

# Display the five first rows 
actual_profit_w_pred_intervals %>% head(5)

# Function to get the number of bets placed, the proportion of winning/losing bets
# and the overall profit for each quantile-based decision rule 

summarise_bets_placed <- function(quantile){
  
  dt <- copy(actual_profit_w_pred_intervals)
  
  # Distinguish positive profits from negative profits  
  dt  <- dt[, positive_profit := ifelse(Y_test > 0, 1, 0)]
  
  dt <-  dt[get(quantile) > 0,
            .(nb_bets_placed = .N,
              sum_profit = sum(Y_test)),
            by = .(positive_profit)]
  
  dt$decision_rule_quantile <- quantile
  
  dt  <- dt[,
            c("prop_bets", 
              "sum_nb_bets_placed",
              "sum_profit_quantile") := list(nb_bets_placed/sum(nb_bets_placed),
                                             sum(nb_bets_placed),
                                             sum(sum_profit)),
            by= .(decision_rule_quantile)]
  
  sdcols <- c("decision_rule_quantile", "positive_profit", "nb_bets_placed",
              "prop_bets", "sum_nb_bets_placed", "sum_profit", "sum_profit_quantile")
  
  dt <- dt[, .SD, .SDcols = sdcols]
  
  return(dt)
  
}

quantiles <- list("Q10", "Q20", "Q30", "Q40", "Q50",
                  "Q60", "Q70", "Q80", "Q90", "Q100")

bets_placed_quantile <- map(quantiles,
                            summarise_bets_placed)

bets_placed_quantile <- rbindlist(bets_placed_quantile)

bets_placed_quantile

# Random sample from test set
set.seed(42) # for reproducibility 

size <- 2500
ind <- sample(1:nrow(X_test), size)
X_test_samp <- X_test[ind]
Y_test_samp <- Y_test[ind]
sum(Y_test_samp)

# Predict on the test set
quantiles_reg_forests <- predict(Model_final, 
                                 data = X_test_samp,
                                 type = "quantiles",
                                 quantiles = seq(0, 1, by = 0.1))

pred_intervals <- as.data.table(quantiles_reg_forests$predictions)

actual_profit_w_pred_intervals <- data.table(cbind(Y_test_samp, pred_intervals))

names(actual_profit_w_pred_intervals) <- c("Y_test","Q0", "Q10","Q20", "Q30", "Q40",
                                           "Q50","Q60", "Q70","Q80", "Q90", "Q100")

quantiles <- list("Q10", "Q20", "Q30", "Q40", "Q50",
                  "Q60", "Q70", "Q80", "Q90", "Q100")

bets_placed_quantile_samp <- map(quantiles,
                                 summarise_bets_placed)

bets_placed_quantile_samp <- rbindlist(bets_placed_quantile_samp)

bets_placed_quantile_samp

# Worst-Case Scenario - we only keep bets with negative profit  
ind <- which(Y_test < 0)
X_test_wc <-  X_test[ind]
Y_test_wc <- Y_test[ind]

# Quuntile Predictions with Random Forest  
quantiles_reg_forests <- predict(Model_final, 
                                 data = X_test_wc,
                                 type = "quantiles",
                                 quantiles = seq(0, 1, by = 0.1))

pred_intervals <- as.data.table(quantiles_reg_forests$predictions)

actual_profit_w_pred_intervals <- data.table(cbind(Y_test_wc, pred_intervals))

names(actual_profit_w_pred_intervals) <- c("Y_test","Q0", "Q10","Q20", "Q30", "Q40",
                                           "Q50","Q60", "Q70","Q80", "Q90", "Q100")

quantiles <- list("Q10", "Q20", "Q30", "Q40", "Q50", 
                  "Q60", "Q70", "Q80", "Q90", "Q100")

bets_placed_quantile_wc <- map(quantiles,
                               summarise_bets_placed)

bets_placed_quantile_wc <- rbindlist(bets_placed_quantile_wc)

bets_placed_quantile_wc